09 - Remote Sensing with Landsat package
Working with satellite imagery in R
In this exercise you learn to use the functions of the
landsat package to process satellite imagery, specifically,
how to:
- streamline histograms so you can combine images
- manipulate multi-band imagery,
- extract and create new data from images, such as NDVI and SAVI
- classify image values using
kmeansunsupervised algorithm to detect similar areas - segment features in satellite imagery
Task 1: Set up your workspace
Start by installing the required packages first.
rasterVis and RColoBrewer are used for
visualisation of rasters, lattice and
latticeExtra for extra graphical utilities, while
landsat provides the imagery samples and rgl
allows for 3D visualisation.
Task 2: Pre-processing of Landsat datasets
Landsat packages offers sample Landsat satellite imagery decomposed
into single bands, and labelled by the month of capture (25 November
2002). You will load and practice raster analysis with these 300x300
pixel samples. Note that each single band image shows as panchromatic
(greyscale) with values ranging from 0 (white) to 255 (black), with the
interim colors being shades of grey. This is unsigned 8 bit imagery. The
number behind the filename (e.g. nov3) refers to the color
band of the image: 1 = Green, 2 = Blue, 3= Red, 4 = Near-infrared. Once
you combine these color bands and plot with plotRGB(), you
can see the true- or false-color imagery depending on which color you
load into which RGB channel.
The images are in the SpatialGridDataFrame format and
need to be converted to raster format before
manipulation.
## load indvidual band image data from landsat package
`?`(nov)
# load band#3 red channel of the image
data(nov3)
plot(nov3)Task 3: Load elevation data and plot it in 3D (Adela to demo
dem in the landsat package denotes a ‘digital elevation
model’. Once you load it you can convert it to Formal Raster Layer with
raster() function and then continue processing using the
usual raster package functions.
plot3D() is a neat function in the
rasterVis package, which opens a separate windows and plots
elevation values in 3D space if present in the raster. You need to close
the window if you want to update or plot another object.
Plot in 3D with raster
RasterVis package function plot3D() loads a neat 3D
viewer, which opens in a new window If you have the raster
library and don’t see the plot3D() function working, enable
the rglwidget().
I have not been able to get an equivalent 3D image with
terra library
Task 4: Load and explore RGB data components
# let's load data for July and explore it
data("july1")
data("july2")
data("july3")
data("july4")
j1 <- rast(july1) # blue
j2 <- rast(july2) # green
j3 <- rast(july3) # red
j4 <- rast(july4) # near-infrared
## check out the image histogram\t
plot(j1)Task 5: Plot RGB image
library(terra)
# Create a true-color composite (R = j3, G = j2, B = j1)
myRGB <- c(j3, j2, j1) # `c()` combines into a SpatRaster with multiple layers
# Create a false-color composite (CIR: NIR = j4, R = j3, G = j2)
myCIR <- c(j4, j3, j2)
### let's see how the NIR, R, and G bands relate (from lattice)
names(myCIR) <- c("NIR", "Red", "Green")
splom(myCIR, varname.cex = 1) # scatter plot matrixTask 6: Manipulate image rendering by histogram stretching and cloud masking
First, let’s use histogram stretch
Next, let’s try linear stretch
## different stretches - here linear stretch
plotRGB(myRGB, stretch = "lin", main = "True Color RGB")Any idea what the red color represents in myCIR?
Finally, drape one of the images over a 3D model. This bit can be a
bit touchy, relies on the raster library and take a while to get to
work. I tend to run the plot3D() lines alone to generate
the 3D view in a pop-up window rather than rmarkdown output. Everytime
you wish to refresh the view, you must close the pop-up
window.
## finally... rglwidget() # this widget helps get the first view rendered in
## rmarkdown, refreshing however is more tedious
# ?rglwidget()
t <- plot3D(dem, col = rainbow(255)) ## you need to close RGL device manually first and then run this line!
# Save it to a file. This requires pandoc
filename <- tempfile(fileext = ".html")
htmlwidgets::saveWidget(rglwidget(), filename)
browseURL(filename)
plot3D(dem, drape = raster(july4)) ## should drape image j4 over DEM, if problematic try in .R script and watch for a pop-up widget. Ask Adela to demo!We have a lot of clouds in the July image - mask them!
RStoolbox comes with a suite of pre-processing
functions, including cloudMask to identify clouds in
optical satellite imagery:
library(landsat)
j1 <- rast(july1) # blue
j2 <- rast(july2) # green
j3 <- rast(july3) # red
j4 <- rast(july4) # near-infrared
# Create a false-color composite (CIR: NIR = j4, R = j3, G = j2, B = j1)
myCIR <- c(j4, j3, j2, j1)
# Rename the NIR, R, and G bands
names(myCIR) <- c("NIR", "Red", "Green", "Blue")# Calculate cloud index
library(RStoolbox)
library(ggplot2)
cldmsk <- cloudMask(myCIR, blue = 4, tir = 1)
ggR(cldmsk, 2, geom_raster = TRUE) # geom_raster plots faster, 2 indicates the mask# mask by threshold, region-growing around the core cloud pixels
cldmsk_final <- cloudMask(cldmsk, threshold = 0.1, buffer = 5)
## plot cloudmask
library(ggplot2)
ggRGB(myCIR, stretch = "lin") + ggR(cldmsk_final[[1]], ggLayer = TRUE, forceCat = TRUE,
geom_raster = TRUE) + scale_fill_manual(values = c("yellow"), na.value = NA)Task 7: Create new information from satellite imagery: NDVI
Let us calculate the Normalized Difference Vegetation Index (NDVI) and see where the vegetation grows most in our Landsat image: Remember the formula for NDVI is: (NIR - RED) / (NIR + RED)
class : SpatRaster
dimensions : 300, 300, 1 (nrow, ncol, nlyr)
resolution : 30, 30 (x, y)
extent : 390045, 399045, 4482105, 4491105 (xmin, xmax, ymin, ymax)
coord. ref. :
source(s) : memory
name : /data/LANDSAT/MASTER/Gorskiout/G20021125_7_band4.geo.asc
min value : -0.3114754
max value : 0.5664336
# uncomment and run in 3D plot3D(dem, drape=ndvi, zfac = 1.5)
### remove values below zero
ndvi[ndvi <= 0] <- NA
plot(ndvi)Now you can see the areas of the highest reflectance and thus most healthy vegetation in November
Task 8: Create new information from satellite imagery: SAVI
SAVI stands for Soils Adjusted Vegetation Index, and this is another calibrated view of the ground.
### another index SAVI (soil adjusted vegetation index)
ndvi <- (n4 - n3)/(n4 + n3)
savi <- (n4 - n3)/((n4 + n3) * 0.25) # with L=1 -> similar to NDVI
### let´s compare visually
par(mfrow = c(1, 2))
plot(savi, main = "SAVI")
plot(ndvi, main = "NDVI")Task 9: Unsupervised Classification with k-means
We would like to isolate and better see the clusters of growth within our image. We will run kmeans function on a composite image in order to cluster data based on similarity or similar groups!
# first, let's select an image and make it into a brick including ndvi
data(nov2)
data(nov1)
n2 <- rast(nov2)
n1 <- rast(nov1)
ndviclass : SpatRaster
dimensions : 300, 300, 1 (nrow, ncol, nlyr)
resolution : 30, 30 (x, y)
extent : 390045, 399045, 4482105, 4491105 (xmin, xmax, ymin, ymax)
coord. ref. :
source(s) : memory
name : /data/LANDSAT/MASTER/Gorskiout/G20021125_7_band4.geo.asc
min value : -0.3114754
max value : 0.5664336
integer(0)
# create a new composite brick out of the available data
myNewBrick <- c(n4, n3, n2, n1, ndvi)
splom(myNewBrick)
Next, run the kmeans classification. Beware that the kmeans function
does not tolerate NA/INF/NaN and similar values. Our new brick should
not have any but in future classification remember that you need to get
around them, either by exclusion or substitution via mean values.
# Run kmeans classification on the values in your new brick Read on
# Thresholding here:
# https://rspatial.org/raster/rs/3-basicmath.html#vegetation-indices
ICE_df <- as.data.frame(myNewBrick)
set.seed(99)
cluster_ICE <- kmeans(ICE_df, 4) ### kmeans, with 4 clusters
str(cluster_ICE)List of 9
$ cluster : Named int [1:90000] 4 3 4 4 4 1 4 4 3 3 ...
..- attr(*, "names")= chr [1:90000] "1" "2" "3" "4" ...
$ centers : num [1:4, 1:5] 48.2 37.5 79.8 59.2 39.9 ...
..- attr(*, "dimnames")=List of 2
.. ..$ : chr [1:4] "1" "2" "3" "4"
.. ..$ : chr [1:5] "/data/LANDSAT/MASTER/Gorskiout/G20021125_7_band4.geo.asc" "/data/LANDSAT/MASTER/Gorskiout/G20021125_7_band3.geo.asc" "/data/LANDSAT/MASTER/Gorskiout/G20021125_7_band2.geo.asc" "/data/LANDSAT/MASTER/Gorskiout/G20021125_7_band1.geo.asc" ...
$ totss : num 20611557
$ withinss : num [1:4] 1042091 848325 839810 1215809
$ tot.withinss: num 3946035
$ betweenss : num 16665522
$ size : int [1:4] 34968 29608 8009 17415
$ iter : int 3
$ ifault : int 0
- attr(*, "class")= chr "kmeans"
# cluster_ICE <- cluster::clara(ICE_df, 4) ### another option, clara, with 4
# clusters
# convert cluster information into a raster for plotting
clusters <- rast(myNewBrick) ## create an empty raster with same extent than ICE
clusters <- setValues(clusters, cluster_ICE$cluster) # convert cluster values into raster
clustersclass : SpatRaster
dimensions : 300, 300, 5 (nrow, ncol, nlyr)
resolution : 30, 30 (x, y)
extent : 390045, 399045, 4482105, 4491105 (xmin, xmax, ymin, ymax)
coord. ref. :
source(s) : memory
names : /data/~eo.asc, /data/~eo.asc, /data/~eo.asc, /data/~eo.asc, /data/~eo.asc
min values : 1, 1, 1, 1, 1
max values : 4, 4, 4, 4, 4
# uncomment to plot the clusters in 3D over the DEM plot3D(dem, drape=clusters,
# col=c('white', 'green', 'blue', 'yellow'))
# calculate the average spectral signature of 1-4 bands of growth
ICE_mean <- zonal(myNewBrick, clusters, fun = "mean")
ICE_mean # see the values for ndvi (layer) being most distinct X.data.LANDSAT.MASTER.Gorskiout.G20021125_7_band4.geo.asc
1 1
2 2
3 3
4 4
X.data.LANDSAT.MASTER.Gorskiout.G20021125_7_band3.geo.asc
1 1
2 2
3 3
4 4
X.data.LANDSAT.MASTER.Gorskiout.G20021125_7_band2.geo.asc
1 1
2 2
3 3
4 4
X.data.LANDSAT.MASTER.Gorskiout.G20021125_7_band1.geo.asc
1 1
2 2
3 3
4 4
X.data.LANDSAT.MASTER.Gorskiout.G20021125_7_band4.geo.asc.1
1 1
2 2
3 3
4 4
/data/LANDSAT/MASTER/Gorskiout/G20021125_7_band4.geo.asc
1 48.23805
2 37.47494
3 79.83131
4 59.23101
/data/LANDSAT/MASTER/Gorskiout/G20021125_7_band3.geo.asc
1 39.87062
2 33.39692
3 41.67599
4 45.38708
/data/LANDSAT/MASTER/Gorskiout/G20021125_7_band2.geo.asc
1 39.84423
2 35.98109
3 45.57810
4 44.90479
/data/LANDSAT/MASTER/Gorskiout/G20021125_7_band1.geo.asc
1 55.70399
2 53.02459
3 58.27282
4 58.88780
/data/LANDSAT/MASTER/Gorskiout/G20021125_7_band4.geo.asc.1
1 0.09473271
2 0.05560162
3 0.31161705
4 0.13208106
Note that you have aggregated the final cluster raster by using the
focal() function using the mean function on
the four clusters identified by kmeans() as similar.
Task 10: What is the trend in de-/afforestation? - Individual tree crown segmentation
The ITC (Individual Tree Crowns) delineation approach finds
local maxima within imagery that contains subtle color
differences, such as the canopy image provided. The
itcIMG() function designates these maxima as tree tops,
then uses a decision tree method to grow individual crowns around the
local maxima.
The image we use is based on LiDAR (Light Detection and Ranging) in
xyz format. To get the segmentation going you will need the
terra and itcSegment packages
# install.packages('itcSegment')
library(terra)
library(itcSegment)
forest <- rast("../data/imgData.tif")
plot(forest)# Use the itcIMG() function to detect and grow the individual crowns
`?`(itcIMG())
se <- itcIMG(forest, epsg = 32632)
summary(se) ID X Y CA_m2
Min. : 1.00 Min. :676744 Min. :5091658 Min. : 6.48
1st Qu.: 30.25 1st Qu.:676754 1st Qu.:5091661 1st Qu.:34.12
Median : 59.50 Median :676776 Median :5091663 Median :43.74
Mean : 59.50 Mean :676782 Mean :5091663 Mean :45.66
3rd Qu.: 88.75 3rd Qu.:676802 3rd Qu.:5091666 3rd Qu.:56.17
Max. :118.00 Max. :676832 Max. :5091669 Max. :95.58
class : SpatVector
geometry : polygons
dimensions : 118, 4 (geometries, attributes)
extent : 676744.4, 676837.1, 5091659, 5091747 (xmin, xmax, ymin, ymax)
coord. ref. : WGS 84 / UTM zone 32N (EPSG:32632)
names : ID X Y CA_m2
type : <int> <num> <num> <num>
values : 1 6.768e+05 5.092e+06 59.94
2 6.768e+05 5.092e+06 82.62
3 6.768e+05 5.092e+06 54.27
### Let´s overlay the image and the product of segmentation (run both lines)
plot(forest)
plot(se, axes = F, add = T, main = "Lidar image of forest with segmented tree-crowns")Task 11: Visualise the segmentation result in Leaflet
You can probably do all of this yourself, but here is a hint about projecting the SpatVector, just in case:
class : SpatVector
geometry : polygons
dimensions : 118, 4 (geometries, attributes)
extent : 676744.4, 676837.1, 5091659, 5091747 (xmin, xmax, ymin, ymax)
coord. ref. : WGS 84 / UTM zone 32N (EPSG:32632)
names : ID X Y CA_m2
type : <int> <num> <num> <num>
values : 1 6.768e+05 5.092e+06 59.94
2 6.768e+05 5.092e+06 82.62
3 6.768e+05 5.092e+06 54.27
[1] "PROJCRS[\"WGS 84 / UTM zone 32N\",\n BASEGEOGCRS[\"WGS 84\",\n ENSEMBLE[\"World Geodetic System 1984 ensemble\",\n MEMBER[\"World Geodetic System 1984 (Transit)\"],\n MEMBER[\"World Geodetic System 1984 (G730)\"],\n MEMBER[\"World Geodetic System 1984 (G873)\"],\n MEMBER[\"World Geodetic System 1984 (G1150)\"],\n MEMBER[\"World Geodetic System 1984 (G1674)\"],\n MEMBER[\"World Geodetic System 1984 (G1762)\"],\n MEMBER[\"World Geodetic System 1984 (G2139)\"],\n ELLIPSOID[\"WGS 84\",6378137,298.257223563,\n LENGTHUNIT[\"metre\",1]],\n ENSEMBLEACCURACY[2.0]],\n PRIMEM[\"Greenwich\",0,\n ANGLEUNIT[\"degree\",0.0174532925199433]],\n ID[\"EPSG\",4326]],\n CONVERSION[\"UTM zone 32N\",\n METHOD[\"Transverse Mercator\",\n ID[\"EPSG\",9807]],\n PARAMETER[\"Latitude of natural origin\",0,\n ANGLEUNIT[\"degree\",0.0174532925199433],\n ID[\"EPSG\",8801]],\n PARAMETER[\"Longitude of natural origin\",9,\n ANGLEUNIT[\"degree\",0.0174532925199433],\n ID[\"EPSG\",8802]],\n PARAMETER[\"Scale factor at natural origin\",0.9996,\n SCALEUNIT[\"unity\",1],\n ID[\"EPSG\",8805]],\n PARAMETER[\"False easting\",500000,\n LENGTHUNIT[\"metre\",1],\n ID[\"EPSG\",8806]],\n PARAMETER[\"False northing\",0,\n LENGTHUNIT[\"metre\",1],\n ID[\"EPSG\",8807]]],\n CS[Cartesian,2],\n AXIS[\"(E)\",east,\n ORDER[1],\n LENGTHUNIT[\"metre\",1]],\n AXIS[\"(N)\",north,\n ORDER[2],\n LENGTHUNIT[\"metre\",1]],\n USAGE[\n SCOPE[\"Navigation and medium accuracy spatial referencing.\"],\n AREA[\"Between 6°E and 12°E, northern hemisphere between equator and 84°N, onshore and offshore. Algeria. Austria. Cameroon. Denmark. Equatorial Guinea. France. Gabon. Germany. Italy. Libya. Liechtenstein. Monaco. Netherlands. Niger. Nigeria. Norway. Sao Tome and Principe. Svalbard. Sweden. Switzerland. Tunisia. Vatican City State.\"],\n BBOX[0,6,84,12]],\n ID[\"EPSG\",32632]]"
# Project the SpatialVector using terra library
se4326 <- project(se, "EPSG:4326")
# ...then we can combine them with leaflet
library(leaflet)
leaflet() %>%
addTiles() %>%
addRasterImage(forest) %>%
addPolygons(data = se4326, weight = 1, color = "black") # TadaaCan you figure out where this forested landscape is from? And how many tree crowns have you just detected in how big an area? What is the average tree density?
Task 12: a more demanding segmentation example with a bigger LIDAR image!
And because we have it in the data/ folder, try with this larger example.
r <- rast("../data/myDem_subset.tif")
r
plot(r)
se_large <- itcIMG(r, epsg = 25829, ischm = T) # can be a bit slow (2-3mins)
summary(se_large)
plot(r)
plot(se_large, axes = F, add = TRUE)
# Adjust 'th' argument for excessive capture of small growth
se_large5 <- itcIMG(r, epsg = 25829, th = 5, ischm = T) # th - how low should algorithm be looking for canopy
# Plot raster with axes (plotting SpatVectors with terra tends to duplicate
# axes)
plot(r, axes = TRUE, main = "Forest LiDar with tree-crown overlay")
# Add vector without drawing axes again
plot(se_large5, border = "white", axes = FALSE, add = TRUE)Again, what is the area of your raster and how many trees or bushes did you detect? Where do you think you are?
Want to see the result in Leaflet?
# Write a SpatVector to a Shapefile writeVector(se_large5,
# '../data/itcTrees_subset.shp', filetype = 'ESRI Shapefile')
# Load shapefile as SpatVector
rse <- vect("../data/itcTrees_subset/itcTrees_subset.shp")
# Reproject to WGS 84 (EPSG:4326)
rse4326 <- project(rse, "EPSG:4326")
# Control question: where is this landscape from?
r <- rast("../data/myDem_subset.tif")
leaflet() %>%
addTiles() %>%
addProviderTiles("Esri.WorldPhysical") %>%
# addProviderTiles('Esri.WorldImagery') %>%
addRasterImage(r) %>%
addPolygons(data = rse4326, weight = 1, color = "black") # Neat :)Task 13: Classify green areas and find the tree crowns in downtown Aarhus
library(RStoolbox)
# unsupervised classification with 3 classes
uc <- unsuperClass(a, nClasses = 4)
# plot result using ggr (ggplot for rasters)
library(ggplot2)
ggR(uc$map, geom_raster = T, forceCat = T) + scale_fill_manual(values = c("grey",
"darkgreen", "sandybrown", "blue"))Try supervised classification of Aarhus landuse with training data
[1] "PROJCRS[\"ETRS89 / UTM zone 32N\",\n BASEGEOGCRS[\"ETRS89\",\n DATUM[\"IRENET95\",\n ELLIPSOID[\"GRS 1980\",6378137,298.257222101,\n LENGTHUNIT[\"metre\",1,\n ID[\"EPSG\",9001]]]],\n PRIMEM[\"Greenwich\",0,\n ANGLEUNIT[\"degree\",0.0174532925199433,\n ID[\"EPSG\",9122]]]],\n CONVERSION[\"Transverse Mercator\",\n METHOD[\"Transverse Mercator\",\n ID[\"EPSG\",9807]],\n PARAMETER[\"Latitude of natural origin\",0,\n ANGLEUNIT[\"degree\",0.0174532925199433],\n ID[\"EPSG\",8801]],\n PARAMETER[\"Longitude of natural origin\",9,\n ANGLEUNIT[\"degree\",0.0174532925199433],\n ID[\"EPSG\",8802]],\n PARAMETER[\"Scale factor at natural origin\",0.9996,\n SCALEUNIT[\"unity\",1],\n ID[\"EPSG\",8805]],\n PARAMETER[\"False easting\",500000,\n LENGTHUNIT[\"metre\",1],\n ID[\"EPSG\",8806]],\n PARAMETER[\"False northing\",0,\n LENGTHUNIT[\"metre\",1],\n ID[\"EPSG\",8807]]],\n CS[Cartesian,2],\n AXIS[\"easting\",east,\n ORDER[1],\n LENGTHUNIT[\"metre\",1,\n ID[\"EPSG\",9001]]],\n AXIS[\"northing\",north,\n ORDER[2],\n LENGTHUNIT[\"metre\",1,\n ID[\"EPSG\",9001]]]]"
Coordinate Reference System:
User input: PROJCRS["ETRS89 / UTM zone 32N",
BASEGEOGCRS["ETRS89",
DATUM["IRENET95",
ELLIPSOID["GRS 1980",6378137,298.257222101,
LENGTHUNIT["metre",1,
ID["EPSG",9001]]]],
PRIMEM["Greenwich",0,
ANGLEUNIT["degree",0.0174532925199433,
ID["EPSG",9122]]]],
CONVERSION["Transverse Mercator",
METHOD["Transverse Mercator",
ID["EPSG",9807]],
PARAMETER["Latitude of natural origin",0,
ANGLEUNIT["degree",0.0174532925199433],
ID["EPSG",8801]],
PARAMETER["Longitude of natural origin",9,
ANGLEUNIT["degree",0.0174532925199433],
ID["EPSG",8802]],
PARAMETER["Scale factor at natural origin",0.9996,
SCALEUNIT["unity",1],
ID["EPSG",8805]],
PARAMETER["False easting",500000,
LENGTHUNIT["metre",1],
ID["EPSG",8806]],
PARAMETER["False northing",0,
LENGTHUNIT["metre",1],
ID["EPSG",8807]]],
CS[Cartesian,2],
AXIS["easting",east,
ORDER[1],
LENGTHUNIT["metre",1,
ID["EPSG",9001]]],
AXIS["northing",north,
ORDER[2],
LENGTHUNIT["metre",1,
ID["EPSG",9001]]]]
wkt:
PROJCRS["ETRS89 / UTM zone 32N",
BASEGEOGCRS["ETRS89",
DATUM["IRENET95",
ELLIPSOID["GRS 1980",6378137,298.257222101,
LENGTHUNIT["metre",1,
ID["EPSG",9001]]]],
PRIMEM["Greenwich",0,
ANGLEUNIT["degree",0.0174532925199433,
ID["EPSG",9122]]]],
CONVERSION["Transverse Mercator",
METHOD["Transverse Mercator",
ID["EPSG",9807]],
PARAMETER["Latitude of natural origin",0,
ANGLEUNIT["degree",0.0174532925199433],
ID["EPSG",8801]],
PARAMETER["Longitude of natural origin",9,
ANGLEUNIT["degree",0.0174532925199433],
ID["EPSG",8802]],
PARAMETER["Scale factor at natural origin",0.9996,
SCALEUNIT["unity",1],
ID["EPSG",8805]],
PARAMETER["False easting",500000,
LENGTHUNIT["metre",1],
ID["EPSG",8806]],
PARAMETER["False northing",0,
LENGTHUNIT["metre",1],
ID["EPSG",8807]]],
CS[Cartesian,2],
AXIS["easting",east,
ORDER[1],
LENGTHUNIT["metre",1,
ID["EPSG",9001]]],
AXIS["northing",north,
ORDER[2],
LENGTHUNIT["metre",1,
ID["EPSG",9001]]]]
[1] "water" "forest" "" "building" "paved"
train <- train[-7, ]
train$Type <- factor(train$Type, levels = c("building", "forest", "paved", "water"))
# plot input data
ggRGB(a, r = 3, g = 2, b = 1, stretch = "lin") + geom_sf(data = train, aes(fill = Type)) +
scale_fill_manual(values = c("yellow", "darkgreen", "sandybrown", "blue"))# fit random forest (splitting training into 70\% training data, 30\%
# validation data)
sc <- superClass(a, trainData = train, responseCol = "Type", model = "rf", tuneLength = 1,
trainPartition = 0.7)
# print model performance and confusion matrix
sc$modelFit[[1]]
TrainAccuracy TrainKappa method
1 0.8880996 0.8472588 rf
[[2]]
Cross-Validated (5 fold) Confusion Matrix
(entries are average cell counts across resamples)
Reference
Prediction building forest paved water
building 667.8 5.6 66.4 23.6
forest 22.6 326.6 9.6 13.2
paved 83.2 2.6 482.4 3.0
water 14.6 3.0 4.6 523.2
Accuracy (average) : 0.8881
# plotting: convert class IDs to class labels (factorize) and plot
r <- as.factor(sc$map)
levels(r) <- data.frame(ID = 1:4, class_supervised = levels(train$Type))
ggR(r, geom_raster = T, forceCat = T) + scale_fill_manual(values = c("sandybrown",
"darkgreen", "yellow", "blue"))Task 14: Optional: Spectral unmixing
Sometimes you want just a two-color/binary output, showing the breakdown on two main landuses (e.g. water vs urban, or forest vs paved, etc.) RStoolbox offers spectral unmixing by implementing the Multiple Endmember Spectral Mixture Analysis (MESMA) approach for estimating fractions of spectral classes, such as spectra of surfaces or materials, on a sub-pixel scale. You can adapt the workflow below for an image of your choice.
The following workflow shows a simple Spectral Mixture Analysis (SMA) with single endmembers per class, extracted from the lsat example image by cell id:
library(RStoolbox)
library(terra)
# to perform a SMA, use a single endmember per class, row by row:
em <- data.frame(lsat[c(5294, 47916)])
rownames(em) <- c("forest", "water")
# umix the lsat image
probs <- mesma(img = lsat, em = em)
plot(probs)Interactively select raster IDs
library(terra)
# Plot the image (pick one band or RGB to visualize)
plotRGB(a, stretch = "lin") # or RGB with plotRGB() if applicable
# Interactively click on the map (e.g., to pick endmembers)
clicks <- click(a, n = 4, cells = TRUE) # 'n' = how many points to click
# View clicked cell IDs and coordinates
print(clicks)The End
Similar approaches can be used when mapping socio-cultural phenomena in satellite imagery, from mosaicing images of night lights as proxies of economic performance, or detecting phenomena in the landscape such as burial mounds, growing urban sprawl, or tracing the outlines of scanned line drawings. (In the latter two you may need to base the classification on reflectance or edge detection rather than elevation.)
References
https://geoscripting-wur.github.io/AdvancedRasterAnalysis/ https://bleutner.github.io/RStoolbox/ http://rspatial.org/spatial/rst/8-rastermanip.html http://neondataskills.org/R/Image-Raster-Data-In-R/ https://geoscripting-wur.github.io/IntroToRaster/ http://wiki.landscapetoolbox.org/doku.php/remote_sensing_methods:home https://rpubs.com/alobo/vectorOnraster